﻿#include "opencv2/opencv.hpp"
#include "opencv2/imgproc/imgproc_c.h"
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"

#include <iostream>
#include <fstream>

using namespace cv;
using namespace std;

int OtsuAlgThreshold(const Mat);

//3. 在OpenCV安装目录下找到课程对应演示图片(安装目录\sources\samples\data)，首先计算灰度直方图，进一步使用大津算法进行分割，并比较分析分割结果。 

/*Otsu实现思路
(1) 计算0~255各灰阶对应的像素个数，保存至一个数组中，该数组下标是灰度值，保存内容是当前灰度值对应像素数
(2) 计算背景图像的平均灰度、背景图像像素数所占比例
(3) 计算目标图像的平均灰度、目标图像像素数所占比例
(4) 遍历0~255各灰阶，计算并寻找类间方差极大值
*/
int main(void)
{
	string path1 = "E:\\AI\\AIsoftware\\opencv\\sources\\samples\\data\\pic2.png";
	string path2 = "E:\\AI\\AIsoftware\\opencv\\sources\\samples\\data\\pic6.png";

	//图pic2.png
	Mat image1 = imread(path1);
	//imshow("原图", image1);

	cvtColor(image1, image1, CV_RGB2GRAY);
	Mat imageOutput1;

	int thresholdValue = OtsuAlgThreshold(image1);
	threshold(image1, imageOutput1, thresholdValue, 255, CV_THRESH_BINARY);
	imshow("pic2.png大津算法", imageOutput1);

	//Mat imageOtsu;
	//threshold(image, imageOtsu, 0, 255, CV_THRESH_OTSU); //Opencv Otsu算法
	//imshow("Opencv Otsu", imageOtsu);
	waitKey(1000);

	//图pic6.png
	Mat image2 = imread(path2);
	//imshow("原图", image2);

	cvtColor(image2, image2, CV_RGB2GRAY);

	Mat imageOutput2;
	int thresholdValue2 = OtsuAlgThreshold(image2);
	threshold(image2, imageOutput2, thresholdValue2, 255, CV_THRESH_BINARY);
	imshow("pic6.png大津算法", imageOutput2);

	waitKey(0);
}

int OtsuAlgThreshold(const Mat image)
{
	if (image.channels() != 1)
	{
		cout << "请读取灰度图！\n" << endl;
		return 0;
	}
	int T = 0; //Otsu算法阈值
	double varValue = 0; //类间方差中间值保存
	double w0 = 0; //目标像素点数所占比例
	double w1 = 0; //背景像素点数所占比例
	double u0 = 0; //目标平均灰度
	double u1 = 0; //背景平均灰度
	double Histogram[256] = { 0 }; //灰度直方图，下标是灰度值，保存内容是灰度值对应的像素点总数
	uchar *data = image.data;
	double totalNum = image.rows*image.cols; //像素总数
	//计算灰度直方图分布，Histogram数组下标是灰度值，保存内容是灰度值对应像素点数
	for (int i = 0; i < image.rows; i++)   //为表述清晰，并没有把rows和cols单独提出来
	{
		for (int j = 0; j < image.cols; j++)
		{
			Histogram[data[i*image.step + j]]++;
		}
	}
	for (int i = 0; i < 255; i++)
	{
		//每次遍历之前初始化各变量
		w1 = 0;
		u1 = 0;
		w0 = 0;
		u0 = 0;

		//背景各分量值计算
		for (int j = 0; j <= i; j++) //背景部分各值计算
		{
			w1 += Histogram[j];  //背景部分像素点总数
			u1 += j * Histogram[j]; //背景部分像素总灰度和
		}
		if (w1 == 0) //背景部分像素点数为0时退出
		{
			continue;
		}
		u1 = u1 / w1; //背景像素平均灰度
		w1 = w1 / totalNum; // 背景部分像素点数所占比例

		//目标各分量值计算
		for (int k = i + 1; k < 255; k++)
		{
			w0 += Histogram[k];  //目标部分像素点总数
			u0 += k * Histogram[k]; //目标部分像素总灰度和
		}
		if (w0 == 0) //从部分像素点数为0时退出
		{
			break;
		}
		u0 = u0 / w0; //目标像素平均灰度
		w0 = w0 / totalNum; // 目标部分像素点数所占比例

		//类间方差计算
		double varValueI = w0 * w1*(u1 - u0)*(u1 - u0); //当前类间方差计算
		if (varValue < varValueI)
		{
			varValue = varValueI;
			T = i;
		}
	}
	return T;
}